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Abstract 

A double-index atomic partitioning of the molecular first-order density matrix is proposed. Contributions 
diagonal in the atomic indices correspond to atomic density matrices, whereas off-diagonal contributions 
carry information about the bonds. The resulting matrices have good localization properties, in contrast 
to single-index atomic partitioning schemes of the molecular density matrix. It is shown that the electron 
density assigned to individual atoms, when derived from the density matrix partitioning, can be made con- 
sistent with well-known partitions of the electron density over AIM basins, either with sharp or with fuzzy 
boundaries. The method is applied to a test set of about 50 molecules, representative for various types of 
chemical binding. A close correlation is observed between the trace of the bond matrices and the SEDI 
(shared electron density index) bond index. 
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I. INTRODUCTION 



The most common way chemists look at molecules is to consider them as composed of atoms 
held together by chemical bonds. Moreover, the chemical characteristics of atoms and functional 
groups of bonded atoms are highly transferable between different molecules. Although this picture 
predates quantum mechanics, it is so useful for rationalizing and even predicting experimental 
observations that it is still ubiquitous. However, the question remains of how to properly describe 
the atom in the molecule (AIM) in a quantum mechanical way. This question has been addressed 
by many people and various techniques have been developed to describe this elusive concept. 

The techniques used thus far can largely be divided into only a few different categories. In 
one group, one uses the attachment of basis functions to atomic centers to extract the AIM. The 
best known method of this type is the Mulliken population analysis JTJ. The second, and for 
the present paper the most important, group of techniques uses a three-dimensional (3D) splitting 
of space with either sharp boundaries between different AIM (e.g. Bader's Quantum Chemical 
Topology (QCT) [0-31), or with more fuzzy boundaries (e.g. the original Hirshfeld method 
and recent extensions [joT-TlOl, and Mayer's fuzzy atoms QTJ). In the second group of methods, 
one uses the molecular electron density and its properties as the guide for obtaining the AIM. 
However, not all AIM properties can be directly expressed in terms of the electron density. A 
very simple example is the kinetic energy of an AIM, for which the full (nondiagonal) density 
matrix is needed rather than the electron density. For some quantities one even has to go up to the 
second order density matrix. This means that a more fundamental approach to the AIM should be 
based on density matrices [fT2|. In Bader's QCT only the electron density and its derivatives are 
required to arrive at the AIM energy via the AIM virial theorem [2J, but several concerns remain. 
For instance, Cioslowski and Karwowski arrived at the conclusion that arbitrary choices in the 
Lagrangian density can have an important influence on the uniqueness of Bader's AIM lfT3l . 

In the present work, we describe how an AIM density matrix that is consistent with a 3D 
division of the molecular density in AIM's can be obtained. Several such methods have been 
explored previously, for instance, in the work of Alcoba et al [fl4]|. These authors derive a QCT- 
based density matrix in the following way. First, the one density matrix is expressed in terms of an 
orthonormal molecular orbital set through the matrix D = {D icrJcr } where ia denotes a molecular 
orbital i with spin a. Given the positive-definite character of this matrix, it can be factorized as 
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follows: 

k,i 

The Kronecker delta 8^,10 can be very simply rewritten as: 

ha,ia = (k<r\la) = (k<r\l(r) A = Yl ( k(7 \W A \l(j) (2) 

A A 

In other words, the Kronecker delta is written as a sum of atom-condensed overlap integrals 
(ka\la) A where the weight W A (r) acts as an operator to delineate the AIM domain in the 
molecule. QCT and Hirshfeld based methods differ mainly in the choice of the operator W A (r), 
which is either binary as in QCT, or fuzzy as in Hirshfeld and related methods. An AIM density 
matrix D^ j(7 with eigenvalues constrained to the interval [0, 1] can then be obtained as: 

< S , = E ( DV2 )„M MWaVv)(D^),^ (3) 

k,l 

This method, although shown to give interesting results when starting from QCT, has as the draw- 
back of being inconsistent with the underlying AIM method. QCT starts from a strict, binary 
division of space in AIM domains. However, inspection of Eq. ([3]) shows that the electron density 
of the AIM is not confined to the AIM domain but spreads over the entire space. Extension of the 
above to a more fuzzy partitioning of space is straightforward, but any scheme along the lines of 
Eq. will result in orbitals extending far outside the atomic basin assigned to A. It should be 
added that Alcoba et al. also introduced a different partitioning by distributing only the molecular 
occupation numbers in the density matrix over the different atoms, retaining the molecular natural 
orbitals. This obviously again does not lead to very well localized density matrices although the 
authors reduce this problem by first localizing the molecular natural orbitals lfT5l[T6Tl . 
In the following we pursue a method that satisfies all the following requirements: 

• The AIM density matrix should be derived starting from a 3D partitioning of space into 
atomic domains such that the AIM density matrix and electron density always remain mu- 
tually consistent. 

• The sum of AIM density matrix eigenvalues should equal the electron occupancy of the AIM 
as obtained from the density analysis, and both the starting AIM density and that obtained 
from the coordinate space diagonal elements of the density matrix should be the same. 

• The density matrices obtained should be localized. 
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• The scheme should involve a double atomic index partitioning with diagonal elements AA 
corresponding to twice the atom A and off-diagonal elements AB, whose eigenvectors cor- 
respond to chemical bonding between the atoms A and B. 

The two-index approach is necessary because of the inherent non-local nature of the density ma- 
trix, as was recently also argued by Mayer and Salvador IfTTII . Introducing the two-index parti- 
tioning also provides an orbital perspective on the changes in the atoms when bonds are formed, 
by extracting bond orbitals with associated eigenvectors from the "bond" density matrices. In this 
sense, our method is reminiscent in philosophy to the so-called Natural Orbitals for Chemical Va- 
lence introduced by Nalewajski et al. lfT8l and used recently to describe chemical bonds by Ziegler 
and co-workers lfT9ll20l 



II. THEORY 



A. Double atom partitioning of the molecular density matrix 

We use notation x = ra to specify the single-electron states in coordinate space, where a 
represents the spin degrees of freedom. The first-order density matrix (1DM) for an iV-electron 
molecule with wave function ^(xx, . . . , x N ) is defined as 



We restrict ourselves to molecules with a singlet ground state. In that case 



(4) 



P(x,x') = -8 a yp(r,r'), (5) 

and the electron spin can be discarded. 

Successful partitioning schemes based on dividing the molecular electron density p(r) = 
p(r,r) in atomic parts rely on introducing positive atomic weight functions W^ifY), which set 
up an atomic decomposition of space [|2TT| in a sense that ^ A Wa{t) = 1. Then, 

p A (r) = p(r)W A (r) (6) 

represents the fraction of the molecular electron density assigned to atom A. The various schemes 
differ mainly in the nature of Wa(v). In Bader's QCT [|2]-[4]|, Wa(^) is a binary operator. In the 
Hirshfeld method [5] and the more recent Hirshfeld-I extension of it Wa(v) is given as: 
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where p A (r) is the density of the isolated atom A. The difference between the regular Hirshfeld 
and the Hirshfeld-I method lies in the choice of the states and charges of the atoms used, as 
described in Bultinck et al. [5]. In the Iterated Stockholder Atoms flU [Kl, a separate weight 
function is used for different spherical shells around the atom. 

In this paper we want to extend the idea of dividing the molecular electron density to the (more 
complex) first-order density matrix. At first, a single atom partitioning scheme 

p(r,r') = J2pA(r,r') (8) 

A 

was sought that would fulfill the essential property of hermiticity with eigenvalues between and 
2. This single atom density matrix should also be strictly localized; that is, the 1DM assigned to 
an atom A should be only appreciably different from zero when both r and r' are near atom A. 
Equivalently, the occupied natural orbitals of the atomic density matrix should be strictly confined 
to the neighborhood of the atom. The transferability of the AIM scheme clearly benefits from this 
concept of localization. Unfortunately, a single atom partitioning scheme of the 1DM as in Eq. ([8]) 
can never fulfill this requirement because the molecular 1DM is (almost always) delocalized and 
has sizeable contributions for r and r' near two different atoms. To accomodate the localization 
requirement there should be a double atomic weighting, depending on both r and r'. This was 
also recently argued by Mayer and Salvador ifTTII . Therefore it seems more natural to introduce a 
double-atom partitioning. 

p(r,r') = Y,PAB(r,r'), (9) 

AB 

where 

Pab{t, r') = - (w A (r)w B (r') + w B (r)w A (r')) p(r, r') = p B A(r, r') (10) 
and the w A (r) are atomic weight functions obeying 

< w A (r) < 1 and ^w A (r) = 1. (11) 

A 

We will use the lowercase notation w A (r) to indicate weight functions used in the double-index 
partitioning ([9 10 ). In general, these will differ from the (uppercase) weight functions W A (r) used 



in a single-index electron density partitioning as in Eq. (|6]). 

Provided the weight functions are properly localized, it is clear that pAB(f, r') also will have 
suitable localization properties, i.e. being appreciably different from zero only when one of the 
arguments is near atom A and the other one near atom B. 



B. Properties of the matrix partitioning 



The individual contributions to the decomposition Q are all hermitian matrices. They clearly 
come in two types, diagonal (A A) and off-diagonal (AB with A=^ B), which will be called atomic 
density matrices and bond matrices, respectively. 

The diagonal terms p AA indeed qualify as first order density matrices, having eigenvalues be- 
tween and 2. The lower bound is a trivial consequence of 

pAA(r, r') = w A (r)p(r, r')w A (r') (12) 

and of the positivity of the molecular 1DM. The upper bound follows from the fact that for any 
electron wave function <f>(r) one has 

J dr J dr , 0(r)p^(r,r')0(r / ) = J dr J dr'[w A {r)(j){r)]p{r, r')[w A {r')<t>{r')] (13) 

< 2 f drw A {rf(j){r) 2 (14) 



< 2 J dr<p{r) 2 . (15) 



The inequality ( fT4| ) results from the upper bound for the eigenvalues of the molecular 1DM, which 
implies that J dr J dr'x(r)x(r')p(r, r') < 2 for any normalized wave function x( r )> an d by ap- 
plying this property to x( r ) = [ w A(f)4>( r )}/{J dr\w A {r)(f){r)) 2 } 1 ^ 2 . The inequality (15) follows 
from w A (r) < 1. 

The summed atomic density matrices do not carry the total number of electrons, since 

^2 / drp AA (r,r) = ^2 drp{r)w 2 A {r) < ^ / drp{r)w A {r) < N; (16) 

A A A 

here we used Eq. (11) and w A (r) < w A (r). The defect in the electron number must be in the bond 
matrices, as Eq. (|9| implies that N = J2ab I drp A B{r, r). For each atom pair AB, the bond 
matrix is seen to carry a positive number of electrons, J drp AB (r, r) = J dr p{r)w A (r)w #(r) > 
0. However, the bond matrices are not positive by construction, and in general negative eigenvalues 
do occur. 

The presence of negative eigenvalues in the bond matrices, disturbing at first, can be readily 



understood by rewriting Eq. ( 10) as 



p AB (r, r') = X><&(rty<&(r') - ^(r)^(r')], (17) 



where ip^Bii 7 ") = \/di/2[w ^(r) ± WB{r))ipi(r), and the ipi(r) and di are the natural orbitals and 
corresponding occupancies of the molecular 1DM. The summed atomic density matrices can be 
rewritten similarly as 

Paa(t, r') + p BB ( ri r') = + ^1(^1^')], (18) 

% 

To see what is going on, consider the most naive picture of covalent binding, with a fully oc- 
cupied (di = 2) molecular orbital ipi(r) m [(f) a^) + 4>b{t)]/\/2 as the bonding combination 
of atomic orbitals 4>a(v) and (f>B( r )- Assuming extreme localization properties for the weight 
functions (i.e. WA{r)(f> B {r) = , WA(r)(f>A(r) = 4>A( r )) one can men interpret the correspond- 
ing ■0. (r) ~ [4>A( r ) =t 0B( r )]/v^2 as the bonding/antibonding combination. In this idealized 



situation, the summed atomic density matrices in Eq. (18) have equal occupancy in the bond 



ing/antibonding orbitals, and the bond matrix in Eq. (17) needs to have negative eigenvalues in 
order to destroy the occupancy of the antibonding combination and enforce the occupation of the 
bonding combination in the total molecular 1DM. Note that negative eigenvalues for AB combi- 
nations also occur in the Natural Orbitals for Chemical Valence (NOCV) technique used recen tly 
by Ziegler and co-workers lfT9ll20l . 



C. Consistency of density matrix and electron-density partitioning 

It is interesting to reflect on what would be the total electron density Pa(^) assigned to a par- 
ticular atom A in the double atomic partitioning scheme of Eq. (|9]). Apart from the diagonal A A 
density matrix, there are now also contributions from the AB bond matrices, and these can be 
distributed over the single-atom densities in different ways. We will study two of these, as they 
can be considered to be extreme cases. 

In the nonweighted scheme (this could also be called a Mulliken-like method Il22l0 . the atoms 
A and B receive an equal share of the density in the bond, 

PaM = PAA(r,r) + pAB(r,r) = p(r)w A {r). (19) 

B&A) 

In the weighted scheme, the atoms receive a share reflecting the balance of the weights wa{v) and 
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w B (r), 



Pa(t) = pAA{r,r)+ pAB{r,r)( 



2w A (r) 



B&A) 



w A (r) + w B (r) 



p(r)w A (r) 



E 

. B 



2w A (r)w B (r) 
w A (r) + w B (r) 



(20) 



Summing Eq. ( |20| ) over all atoms yields the molecular electron density, as it should be: 

p(r) 



5>5(r) 

A 



^ I 2w\(r)w B {r) 



AB 



w A (r) +w B (r] 



p{r) 
p(r) 



1 



^ | 2w A (r)w B (r) + 2w|(r)u; A (r) 



AB 



E 



w A {r) + w B [r) w A (r)+w B (r) 



w A (r)w B (r) jw A (r) + w B (r)) 
w A (r) +w B (r) 



w B {r 



AB 



p{r) 



(21) 



In both schemes, there is now the possibility of choosing the weights w A (r) in such a way that 
the single-atom density p A w (r) coincides with an established AIM electron-density model of the 
form In the nonweighted scheme, one simply takes w A (r) = W A (r). In the weighted scheme, 
consistency requires that the weights obey a set of nonlinear equations, 



WA : w A (r) 



E 

, B 



2w A (r)w B (r) 
w A (r) + w B (r) 



W A (r) 



(22) 



In practice, we always found that the iterative sequence 



1/2 



wf(r) = W A (r); w^'(r) 



(i+l), 



W A (r) 



E 



(23) 



converges rapidly (in at most 20 iterations) to a stable solution of Eq. (22) with 



E 



(24) 



as the convergence criterion. 

It should be noted that in all numerical work in this paper the Hirshfeld-I weight functions 
W A (r) were taken as input; for other "fuzzy atom" prescriptions, like ISA, we expect similar 
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results. The Bader AIM concept is fundamentally different in that it has hard boundaries for the 
atomic basins. In the absence of non-nuclear attractors space is partitioned into the atomic basins, 
and the corresponding Bader weight function Wa(v) is zero/one if r is outside/inside the basin of 
atom A. One can easily verify that for such binary weight functions the nonweighted and weighted 



schemes coincide, and there is no need to solve the nonlinear equations (22). On the other hand, 
the possible presence of non-nuclear attractors is problematic, as it cannot easily be reconciled 
with the underlying AIM picture. However, it deserves to be mentioned that with the present 
density matrix partitioning, the atomic electron densities from QCT remain localized to the AIM 
basin. 



D. Expressions in a finite basis set 

In practice, the partitioned density matrices are expressed in a finite basis set, as used in the 
molecular calculation. In the basis of the molecular natural orbitals ipi(r), e.g., the partitioned 
1DM becomes 

{pAB)ij = J $rd i r''ip i {r)p A B{r,r')'ilj j {r') 

3 rd 3 r'4> i (r)4> j (r')^lwA(r)w B (r') + w B {r)w A {r')} d k i) k (r)i) k (r') 

k 

= ^ d k -[C^Cf k + Cf k Cf k \. (25) 

k 

Note that in the present calculation the tpi(r) are just the molecular HF orbitals, and the natural 
occupations di are 2 (0) for the occupied (unoccupied) orbitals. It follows that the pab matrix is 
expressed in terms of the atomic overlap matrix elements (AOM) : 

C# = J d 3 r^(r> A (r)^-(r). (26) 

The traces of the atomic and bond density matrices correspond exactly to the net and overlap 
populations for fuzzy atoms as defined in Eq. (6) of Ref. [fTTI (of course, when the same weight 
functions are used). 



The matrix representation (pAB)ij in Eq. (25), in the finite single-particle space spanned by 



the basis set, is only equivalent to Eq. ( |T0| ) in the limit of a complete basis set. But since the 
underlying SCF calculation is performed in the same finite basis set, it can be argued that consistent 



calculations actually require the use of the limited basis set expressions in Eq. (25). In a sense, 



the single-particle space spanned by the finite basis set is "all there is". This holds in particular 
for subsequent AIM energy considerations: care should be taken that only matrix elements of 
the electronic Hamiltonian in the finite basis set are used, since these are the ones that fixed the 
electronic structure of the molecule. 

For a reasonable large basis set, the differences between results obtained using Eq. d25j) and 



Eq. (10) are small anyway, as can be suspected by the clean convergence behaviour in Table III 



HI. COMPUTATIONAL METHODS 

The partitioning scheme described in Sec. [TT] was tested by partitioning the 1DM of a small 
set (listed in Table [I]) of ca. 50 simple molecules with a singlet ground state, representative of a 
diverse variety of chemical bonds. The 1DM was calculated at the Hartree-Fock level of theory 
using the Aug-cc-pVDZ basis set [|23] - |25l . The geometry was taken from a B3LYP Il26l - |29l Icc- 
pVDZ fl2ll25l|30l optimization. 

The scheme was implemented using atomic weight functions Wa(t) from a Hirshfeld-I anal- 
ysis. The weight functions ^(r) for the double atomic partitioning in Eq. ^ were constructed 



in both the nonweighted scheme (u>yi(r) = Wa(v)) and in the weighted scheme of Eq. (22). The 
iterative Hirshfeld weights and the adapted weights ^(r) are calculated on atom-centered grids, 
using a logarithmic radial grid of 100 points with r min = 10~ 6 A and r max = 20 A , and Lebedev 
angular grids [f3TT - [36]| with 170 points having randomized orientation on the different shells. 

For CO and C3H3N additional calculations were performed in larger (Aug-cc-pVTZ and Aug- 
cc-pVQZ ) ll23l - l25l basis sets and with larger integration grids, in order to assess basis set and grid 
convergence. 



IV. RESULTS AND DISCUSSION 
A. Natural orbitals and populations 

Diagonalisation of the matrices in Eq. ( [25] ) leads to the natural orbitals and occupations. As 
an example of the results obtained, Figs. [T]{2] depict the dominant natural orbitals and occupations 
of the atomic density matrix paa of carbon and oxygen in CO, within the weighted scheme. The 
natural orbitals are slightly deformed versions of the typical atomic Is (a), 2s (b), 2p x , 2p y and 
2p z (c-d-e) orbitals (with the z-direction corresponding to the internuclear axis). They are clearly 
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CeH6 (benzene) 


CHONH2 


H 3 0+ 


LiF 


OS 2 


A1H 3 


CF4 


CHOOH 


HC1 


LiH 


SF 6 


B2Hg 


CH 2 NH 


Cl 2 


HCN 


LiOH 




BeH 2 


CH 2 


CO 


HCOOCH3 


N 2 




BH 3 


CH2O2 (dioxirane) 


co 2 


HF 


N 2 




C2H2 (acetylene) 


C 3 H 8 (propane) 


F 2 


HN0 2 


NaCl 




C2H4 


CH3NH2 


H 2 


HNO3 


NaOH 




C2H6 


CH3OCH3 


H 2 


HOC1 


NH 3 




C3H3N (acrylonitrile) 


CH3OH 


H 2 S 


HOOH 


2 




C3H4 (cyclopropene) 


CH 4 


H2SO4 


Li 2 


3 





TABLE I: List of molecules in the test set. 



localized and have populations between and 2. Note that for these main contributions the shape 
of the orbitals in the nonweighted scheme is visually indistinguishable from those presented in 
Figs. [T]j2j The orbitals that represent the core Is atomic orbitals have a population well below 
2.00 (the expected value of an orbital not involved in bonding), resp. 1.84 in carbon and 1.90 
in oxygen. The rather low eigenvalues for the essentially core orbitals are less pleasing from a 
chemical point of view. The occupations in the nonweighted scheme are quite different, and are 
indicated between brackets in Figs. [TJj2j In the nonweighted scheme, the populations are more in 
accordance with the full occupation expected of core orbitals, i.e. 1.96 and 1.99 for the carbon 
and oxygen Is orbitals. As discussed later, this difference between both schemes is typical, as 
the weighted scheme tends to assign a larger fraction of electrons to the bond matrices than the 
nonweighted scheme. Hence the nonweighted scheme has invariably larger populations for the 
atomic density matrix natural orbitals. The 2s orbitals have shifted away from the bonding region, 
to accomodate the free electron pair on C and O in the CO molecule. However, they are not fully 
occupied. The 2p x , 2p y and 2p z orbitals (not expected to be fully occupied) are shifted to the 
bonding region. 
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d = 1.842 (1.963) d = 1.579 (1.821) d=0.397 (0.420) d=0.343 (0.372) d=0.343 (0.372) 




FIG. 1: The dominant natural orbitals, and the corresponding occupations, of the atomic density matrix of 
carbon in a CO molecule calculated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The 
occupations in the nonweighted scheme are also indicated, between brackets. 



d = 1.904 (1.988) d = 1.677 (1.894) d = 1.105 (1.276) d = 1.105 (1.276) d=0.905 (1.051) 




FIG. 2: The dominant natural orbitals, and the corresponding occupations, of the atomic density matrix of 
oxygen in a CO molecule calculated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The 
occupations in the nonweighted scheme are also indicated, between brackets. 

Figure [3] depicts the naturals and occupations of the CO bond matrix pab in the CO molecule. 
Only the main contributions are shown. Both schemes result in natural orbital shapes that are very 
similar. There are mainly typical bonding orbitals a (a), tt±, 7r 2 (b-c) with positive eigenvalues, and 
antibonding orbitals a* (h), ix\, -k\ (f-g) with negative eigenvalues. The negative eigenvalues delete 
antibonding contributions of the paa + Pbb matrix, whereas the positive eigenvalues reinforce 
its bonding contributions. The picture also shows two rather "nonbonding" orbitals (d-e) with 
positive eigenvalues. These serve to reinforce the population of the free electron pairs on carbon 
and oxygen (the shifted 2s orbital on these atoms). Note that in the weighted scheme orbitals (i-j) 
appear which correspond to the Is core orbitals on C and O and have a nonnegligible (0. 17 - 0. 10) 
occupation in the bond matrix. 
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d=0.949 (0.929) d^0.892 (0.865) d=0.892 (0.865) d=0.445 (0.270) d=0.345 (0.153) 




FIG. 3: The dominant eigenvectors and occupations of (twice) the CO bond matrix in a CO molecule calcu- 
lated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The occupations in the nonweighted 
scheme are also indicated, between brackets. 

Apart from the dominant contributions shown in Figs.[T]j3j various orbitals with much smaller 
populations are also present (as the Aug-cc-pVDZ molecular basis set has 46 basis functions). A 
complete overview is given in Table |TTJ. Note that for the atomic density matrices, apart from the 
five dominant natural orbitals, only two more have small (~ 1CT 2 — 1CT 3 ), while the remainder 
have vanishing (< 1CT 13 ) populations. For the bond matrix, apart from the ten dominant orbitals, 
only four more have small occupations. This can be understood from the fact that the molecule is 
treated at the HF level, with only 7 (double occupied) spatial orbitals. One can rewrite the atomic 
density matrix as 

N/2 

pAA(r,r') = J~] [w A (r)ipj(r)) [w A {r')ijji(r')] 

i=l 

N/2 

= £^(^(0, (27) 

ii'=i 

where T A is the overlapmatrix of the nonorthogonal basis functions [uu(f )^M r )] 

T$ = J w A {r)A{r)^) d 3 r (28) 
and the ipj(r) form a set of orthonormal basis functions 

J 3 {r) = J2( TA '^ [Mr)A{r)\ ■ (29) 
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C,C 0,0 C,0 

label eigenvalue eigenvalue eigenvalue 

(a) 1.842 ( 1.963 ) 1.904 ( 1.988 ) 0.474 ( 0.464 ) 

(b) 1.579 ( 1.821 ) 1.677 ( 1.894 ) 0.446 ( 0.432 ) 

(c) 0.397 ( 0.420 ) 1.105 ( 1.276 ) 0.446 ( 0.432 ) 

(d) 0.343 ( 0.372 ) 1.105 ( 1.276 ) 0.223 ( 0.135 ) 

(e) 0.343 ( 0.372 ) 0.905 ( 1.051 ) 0.173 ( 0.076 ) 

(f) 0.023 ( 0.008 ) 0.042 ( 0.021 ) -0.170 ( -0.256 ) 

(g) 0.002 ( 0.000 ) 0.005 ( 0.001 ) -0.170 ( -0.256 ) 

(h) < 10~ 13 < 10~ 13 -0.149 ( -0.221 ) 

(i) : : 0.084(0.024) 
(j) 0.050 ( 0.007 ) 

: -0.034 ( -0.061 ) 

-0.007 ( -0.007 ) 
-0.001 ( -0.002 ) 
-0.000 ( -0.000 ) 
< 10~ 13 

Sum 4.528 ( 4.955 ) 6.743 ( 7.508 ) 1.364 ( 0.768 ) 

TABLE II: All natural populations in the atomic density matrices (CC and 00) and bond matrix (CO) in a 
CO molecule calculated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The 
occupations in the nonweighted scheme are indicated between brackets. The labels in the first 
column correspond (for the dominant orbitals) to the labels in Figs.fTpl 



It is clear that pAA( r , r ') is a matrix of rank N/2, since diagonalisation of the N/2 x N/2 overlap- 
matrix T A will yield N/2 eigenvalues. The same reasoning holds for the bond matrix Pab{ t i t'), 



18) with 



but now starting from the 14 (linear independent) functions ipABi( r ) defined in Eq. ( 17 
di = 2forz= l,2,..JV/2. 

The difference between the nonweighted and weighted scheme is again clear from Table |n| The 
differences are sizeable, with consistently smaller populations in the weighted scheme, resulting in 
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traces for the atomic density matrices that are significantly (0.4 - 0.7 ) smaller. This is compensated 
for by larger eigenvalues in the bond matrix in the weighted scheme, for the orbitals not involved 
in bonding. The bonding orbitals (a-b-c) in the bond matrix, however, are about equally populated 
in both schemes. 



CO 


grid=100-170 
Aug-cc-pVDZ 


grid=500-590 
Aug-cc-pVDZ 


grid=100-170 
Aug-cc-pVTZ 


grid=100-170 
Aug-cc-pVQZ 


A , B 


Tt(pa,b) 


Tr(pA,B) 


Tt(pa,b) 




C ,C 
0,0 
CO 


4.528 (4.955) 
6.743 (7.508) 
1.364 (0.768) 


4.528 (4.955) 
6.743 (7.508) 
1.364 (0.768) 


4.529 (4.956) 
6.751 (7.516) 
1.360 (0.764) 


4.531 (4.958) 
6.752 (7.518) 
1.358 (0.762) 



TABLE III: Number of electrons present in the atomic density and bond matrices for CO, in the non- 
weighted and weighted scheme. Bracketed values correspond to the nonweighted scheme. 
The second column contains the results of the default calculation. In the third column the 
number of radial and angular grid points is increased to 500 and 590, respectively. In the 
fourth column a triple-, rather than double-^ basis is used, in the fifth colum a quadruple-^. 



In Table III the stability of the proposed density matrix partitioning is examined. The results 
for the matrix traces are clearly converged as far as grid size is concerned, with deviations of less 
than 0.001 in electron number. The results also seem to be quite stable with respect to basis set 
size, with differences less than 0.01 going from DZ to TZ, and less than 0.002 going from TZ to 
QZ. 




H(4) H(7) 



FIG. 4: Topology of the C3H3N (acrylonitrile) molecule. 



As a second example we analyze some results obtained for the density matrix partitioning in 
C3H3N (acrylonitrile). The topology of the molecule is presented in figure |4| Figure [5] shows the 
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dominant eigenvectors and occupations of the CN bond matrix. The features are very similar to 
those of the CO bond matrix. The bond orbitals o (a), iti (b) and n 2 (c) are about equally populated 
in both schemes. The nonbonding orbital (d) reinforces the population of the free electron pair on 
nitrogen. The antibonding orbitals n\ (e), ti^ (f), and a* (g) have negative eigenvalues that delete 
antibonding contributions of the pc(i)C(i)+PN(2)N{2) matrix. Small but nonnegligible contributions 
exist that are complementary to the Is core electrons on C(l) and N(2) (h and j). All contributions 
that correct the sum of the atomic density matrices (by deleting antibonding parts and reinforcing 
nonbonding parts) are significantly larger within the weighted scheme. 



d=0.899 (0.955) d=0.819 (0.893) d=0.814 (0.881) d=0.335 (0.167) d=-0.301 (-0.493) 




d=-0.290 (-0.479) d=-0.239 (-0.380) d=0.149 (0.048) d=0.129 (0.083) d=0.113 (0.021) 




FIG. 5: The dominant eigenvectors and occupations of the CN bond matrix in a C3H3N molecule calcu- 
lated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The occupations in the nonweighted 
scheme are also indicated, between brackets. 

Figure [6] shows the dominant eigenvectors and occupations of the C(3)-C(5) bond matrix. As 
it is formally a double bond, large bonding orbitals a (a) and 7Ti (b) are expected, next to some 
smaller antibonding contributions a* (d) and 7rJ (e). There are a lot of small contributions not 
shown here involving the 2s orbitals on both carbons, the oc-h bonds and the Is orbitals on both 
carbons, that are complementary to the main orbitals of other atomic and bond matrices. However, 
beyond the many small contributions, there is one considerably larger than the others: a tx 2 (c) 
bond. It's important to mention that this is a rather general feature, also noticed for example in the 
CH bond matrices of CH 4 , where both iri and 7r 2 contributions are important, although there are 
no 7r bonds in the CH 4 molecule. Notice that the trace within the weighted scheme is still larger 
than the one within the nonweighted scheme, while the bonding orbitals of the weighted scheme 
have a significantly lower occupation than that of the nonweighted scheme. 
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FIG. 6: The dominant eigenvectors and occupations of the C(3)-C(5) bond matrix in a C3H3N molecule 
calculated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The occupations in the non- 
weighted scheme are also indicated, between brackets. 

Figure [7] shows the dominant eigenvectors and occupations of the C(l)-C(3) bond matrix. As it 
is formally a single bond, a large bonding orbital a (a) is expected, next to a smaller antibonding 
contribution a* (d). As in the previous case of the C(3)-C(5) double bond, there are a lot of smaller 
contributions, with two of them predominant: a 7Ti (b) and a %2 (c) orbital. Again, the trace within 
the weighted scheme is larger than the trace within the nonweighted scheme, while the bonding 
orbitals have a significantly lower population. 

d=0.725 (0.869) d=0.282 (0.315) d=0.221 (0.229) d=-0.218 (-0.375) 

(a) (b) (c) (d) 

FIG. 7: The dominant eigenvectors and occupations of the C(l)-C(3) bond matrix in a C3H3N molecule 
calculated in the weighted scheme, at the HF/Aug-cc-pVDZ level. The occupations in the non- 
weighted scheme are also indicated, between brackets. 




B. Correlation with shared electron density indices 

A global test for the partitioning scheme is the evaluation of the total population in its bond ma- 
trices. This population should correlate somehow with the bond order. The classical definition of 
bond order equals it to the electronic occupancy of the bonding orbitals minus that of the antibond- 
ing orbitals divided by two. At the Hartree-Fock level of theory, there is a remarkable similarity 
between these classical bond orders and the results of so-called shared electron density indices 
(SEDI) [|TTl[37l - l45l . SEDI are obtained from integration of the exchange density at the Hartree- 
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Fock level or exchange-correlation density at correlated levels of theory. The exchange-correlation 
density reads 

p xc ^ r /) =p (r)p ( r ') - P (2) r '°'\ r '°') ( 30 ) 

where the diagonal elements of the second-order density matrix appear. At the single- 
determinant level, this simplifies to: 

p xc (r,r') = ^p(r,r')p(r',r) (31) 

Integrating r and r' over the atomic domains of atoms A and B and multiplying by two to account 
for the symmetrical integration over B and A eventually results in: 

SEDI(A,B)=AY J Sf J Sf l (32) 

where i and j are occupied molecular orbitals, and S£j and are the atom condensed overlaps. In 
the case that the atomic domains are defined by Hirshfeld-I, the atom condensed overlaps reduce 
to an expression formally very similar to Eq. ([26]). 



Sf j = J dh^{r)W A {r)^{r) (33) 

The main advantage of SEDI is that these can be computed also between non-covalently bonded 
atoms, where the classical expression can no longer be used. Given the classical expression for 
bond order that uses occupancies of bonding and antibonding orbitals, an intriguing question is 
whether the trace of the bond matrices would also give similar results. In order to answer this 
question not only for covalently bonded atoms but in general, this section examines the correlation 
between the trace of bond matrices and the SEDI, both computed starting from the Hirshfeld-I 



analysis. For the SEDI, the Sfi were constructed using the Hirshfeld-I weights Wa(t) in Eq. (33 ). 



The Hirshfeld-I weights were also used to calculate the atomic overlap matrices Cf, of Eq.(26) for 



the nonweighted scheme. The weights Wa{t) of the nonlinear equations of (22 ) were used to build 
the atomic overlap matrices for the weighted scheme. 
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TABLE IV: The SEDI index (second column) versus the traces of the pab matrices in the weighted scheme 
(third column) and the nonweighted scheme (fourth column) for the C3H3N (acrylonitrile) 
molecule. AB contributions corresponding to bonded atoms in the Lewis structure are in bold- 
face. 



Table IV shows a complete comparison between the SEDI index and the matrix traces within 
both the weighted and the nonweighted scheme for the C3H3N (acrylonitrile) molecule. It appears 
that for the triple, double and single bonds present in the molecule (indicated in boldface), the 
SEDI index is quite close to the classical bond order. The matrix traces of the weighted scheme 
are usually somewhat lower, while the matrix traces of the nonweighted scheme are much lower 
(in some cases more than 45 percent). For pairs of atoms that have no classical bond between them, 
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the matrix traces of the weighted scheme are noneglegible and generally somewhat larger than the 
SEDI indices, whereas the matrix traces of the nonweighted scheme are quite small (< 0.08). The 
largest SEDI indices and matrix traces for nonbonded atom pairs can be found in allylic places. 

For the test set of table [j] it was investigated how well the nondiagonal Shared Electron Distri- 
bution Index (SEDI) correlates with the matrix traces within both the weighted and nonweighted 
scheme. A correlation conceptually analogous was first given by I. Mayer for the Mulliken case 
ll46ll . The correlation plots are shown in Fig. [8] For both schemes there is a strong linear correlation 
(R 2 > 0.96), but the proportionality factor is surprisingly close to unity for the weighted scheme 
(0.97). By comparison, the slope for the nonweighted scheme is much smaller (0.60). In the case 
of Hirshfeld-I and the weighted scheme, the SEDI index and the "overlap population" seem to 
give largely equivalent values SEDI(A, B) ~ 2Tr(p^ B ). In view of the extent (50 molecules) 
and diversity of the test set, this can hardly be coincidental. One may argue that this is somewhat 
skewed because of the large number of classically nonbonded atom pairs, for which both measures 
are close to zero. However, restricting to the bonded atom pairs the proportionality factors hardly 
change (0.95 for the weighted, 0.59 for the nonweighted scheme). This is all the more surprising, 
as one would intuitively expect a closer relation between SEDI(A, B) and 2Tr(p AB ) since the 
weight functions are the same (i.e. Wa{t = Wa(v)) for the nonweighted scheme. 




1 1.5 
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FIG. 8: Correlation between the SEDI index and the traces of the bond matrices for the weighted scheme 
(a) and the nonweighted scheme (b) 



The fact that the slope of the correlation plot is larger in the weighted scheme, can be traced 
back to the fact that the bond matrices pab nave a larger trace in the weighted than in the non- 



weighted scheme, as is evident from Table IV The opposite holds for the trace of the atomic den- 
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sity matrices paa' these are larger in the nonweighted scheme than in the weighted scheme. Both 
observations can be explained by the Hirshfeld-I weight functions Wa{t) to be more strongly lo- 



calized around atom A than the weight function ^(r) resulting from the nonlinear equations (22 ). 



This is indeed what is found numerically, and it can also be understood by analyzing Eq. (22): in 
the vicinity of atom A, uu(r) dominates over the weight functions wb{t) of all other atoms. 



Replacing in the denominator of Eq. (22) by wa(v) > wb{t) therefore leads to 

" V Mr) + Mr) V Mr) + Mr) ~ a[ Y ( } 

So the Hirshfeld-I weights Wa(v) are more localized on the individual atoms and will lead to a 
smaller number of shared electrons. 

From the preceding discussion it's clear that there seems to be some freedom in choosing an 
appropriate scheme with corresponding atom weights: one can increase the number of electrons 
in the overlapmatrix pab r ') and get it even close to the SEDI index, at the risk of including 
some parts (e.g. the core Is electrons) that do not belong there. 



V. CONCLUSIONS 



We have introduced a succesful approach to partitioning of a molecular density matrix in con- 
stituent atomic and bond contributions. The partitioning is such that the density matrices and 
electron densities for the atoms in the molecule are mutually consistent, and follows from requir- 
ing that the atom and bond orbitals, as eigenfunctions of the corresponding density matrices, are 
strictly localized. This prevents the use of single atom density matrix partitioning as the molecular 
density matrix is inherently delocalized. 

The atomic density matrices correspond to diagonal elements pAA( r , r ), whereas bond matrices 
correspond to the off-diagonal contributions pab(^, f). Only in the case of the diagonal elements 
are the eigenvalues restricted to the interval [0,2]. For the bond matrices, negative eigenvalues can 
and do occur. 

The weight functions used for partitioning the molecular density matrix were constructed in 
two different ways: using either directly the weights produced from a regular atoms-in-molecules 
(AIM) density based theory (here Hirshfeld-I), or using a weighted scheme. In both cases the 
AIM density derived from the density matrix is equal to the one obtained directly from the AIM 
algorithm. 
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The trace of the bond density matrices was suggested as a useful source of bond indices, loosely 
related to classical bond orders. A remarkably good correlation with shared electron density in- 
dices (SEDI) was found, establishing the chemical relevance of the density matrix partitioning. 
The traces of the bond density matrices can in this way be used quite effectively to characterize 
chemical bonds without requiring second-order density matrices. 

It should be mentioned that while the present analysis is restricted to spin singlet molecules, 
an extension to higher-spin (S > 0) molecular states seems entirely possible, be applying the 



decomposition in Eq. ( 10) to the spin up and down electron density separately, 



p AB (r, r') = l - {w A {r)w B {r') + w B {r)w A {r')) p°{r, r') = p BA (r, r') (35) 

where the atomic weight functions are taken as spin independent. Certainly for Hirshfeld-like 
schemes this seems the most natural choice, as the spin-state of the isolated atoms building up the 
weight functions is washed away when these are considered as molecular constituents. Note that 
even for the simpler problem of molecular electron density partitioning, spin considerations have 
hardly been studied and should be explored further. 
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